

def find_changes(seq, ref):
    """
    Finds all positions in a sequence when the value changes
    from ref to a different value and vice versa.
    """
    changes_to = []
    changes_from = []
    for i in range(1, len(seq)):
        if seq[i - 1] != seq[i]:
            if seq[i] == ref:
                changes_to.append(i)
            if seq[i - 1] == ref:
                changes_from.append(i)
    return changes_to, changes_from


def similarity(seq_a, seq_b):
    """
    Returns percentage of equal values between seq_a and seq_b
    """
    if len(seq_a) != len(seq_b):
        raise ValueError
    count = 0
    for a, b in zip(seq_a, seq_b):
        if a != b:
            count += 1
    return 1 - count / float(len(seq_a))


def closest_left(seq, s, ref):
    c = s // 2
    if c != s and seq[c] > ref:
        return closest_left(seq, c, ref)
    for i in range(s, -1, -1):
        if seq[i] <= ref:
            return i
    return None


def closest(seq, ref):
    """
    Returns closest value to ref in vector seq
    """
    l = closest_left(seq, len(seq) - 1, ref)
    if l is None:
        return 0
    r = min(l + 1, len(seq) - 1)
    return l if abs(ref - seq[l]) <= abs(seq[r] - ref) else r


def hypo(seq, nbins=50):
    from scipy import stats
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib.patches as patches
    import matplotlib.path as path

    print(['z-score = %6.4f' % i for i in stats.zmap([-1, 0, 1], seq)])

    fig = plt.figure()
    ax = fig.add_subplot(111)

    n, bins = np.histogram(seq, nbins)
    left = np.array(bins[:-1])
    right = np.array(bins[1:])
    bottom = np.zeros(len(left))
    top = bottom + n

    # we need a (numrects x numsides x 2) numpy array for the path helper
    # function to build a compound path
    XY = np.array([[left, left, right, right], [bottom, top, top, bottom]]).T

    # get the Path object
    barpath = path.Path.make_compound_path_from_polys(XY)

    # make a patch out of it
    patch = patches.PathPatch(barpath, facecolor='blue', edgecolor='gray', alpha=0.8)
    ax.add_patch(patch)

    # update the view limits
    ax.set_xlim(left[0], right[-1])
    ax.set_ylim(bottom.min(), top.max())
    ax.set_xlabel('distancia')
    ax.set_ylabel('quantitat')

    plt.show()


if __name__ == '__main__':
    averages = []
    for i in range(1, 61):
        with open('dc{0}.states'.format(i)) as fa:
            fa.readline()
            seq_a = [int(x) for x in fa.readline().split()]

        with open('dc{0}.testvit'.format(i)) as fb:
            fb.readline()
            seq_b = [int(x) for x in fb.readline().split()]

        ref = seq_a[0]

        count_a = sum([x - 1 for x in seq_a])
        count_b = sum([x - 1 for x in seq_b])

        print("A: {0} ones {1} twos".format(len(seq_a) - count_a, count_a))
        print("B: {0} ones {1} twos".format(len(seq_b) - count_b, count_b))

        a_to, a_from = find_changes(seq_a, ref)
        b_to, b_from = find_changes(seq_b, ref)

        print("A: {0} changes".format(len(a_to) + len(a_from)))
        print("B: {0} changes".format(len(b_to) + len(b_from)))

        print("similarity: {0} %".format(similarity(seq_a, seq_b) * 100))

        b_to_offset = []
        for b in b_to:
            i = closest(a_to, b)
            b_to_offset.append(b - a_to[i])

        r1 = sum(b_to_offset) / float(len(b_to))

        b_from_offset = []
        for b in b_from:
            i = closest(a_from, b)
            b_from_offset.append(b - a_from[i])

        r2 = sum(b_from_offset) / float(len(b_from))

        changes = b_from_offset + b_to_offset

        avg = sum(changes) / float(len(changes))

        averages.append(avg)

    hypo(averages, 5)
